#!/usr/bin/perl
# 2014/01/17 Batch treatment of BWA, samtool

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use File::Basename;
use Carp qw(croak);
$| = 0;


if($ARGV[0] eq '--help'){
	my $help =<<"EOF";
Usage : $0 -i [sanger fastq file] -j [jump position] -o [old illumina file]
EOF
	die "$help";
}

my $FILE_reference = "/home/hideki/bio/bwa/pombe.ASM294v2.19";	# reference for BWA

my %opt;
getopts("i:j:o", \%opt);
my $FILE_fastq = $opt{i};
my ($name, $dir, $ext) = &fileparse($FILE_fastq, '\..*');
(my $FILE_sai = $FILE_fastq) =~ s/\.fastq/.sai/;
(my $FILE_sam = $FILE_fastq) =~ s/\.fastq/.sam/;
(my $FILE_bam = $FILE_fastq) =~ s/\.fastq/.bam/;
(my $FILE_alnLog = $FILE_fastq) =~ s/\.fastq/_aln.log/;
(my $FILE_samLog = $FILE_fastq) =~ s/\.fastq/_sam.log/;


&status("Start BWA");
#---------------------------------------
# aln
#---------------------------------------
unless(exists $opt{o}){
	system("bwa aln -t 2 $FILE_reference $FILE_fastq > $FILE_sai 2> $FILE_alnLog") == 0 or die "aln command failed: $!";
}else{
	system("bwa aln -t 2 -I $FILE_reference $FILE_fastq > $FILE_sai 2> $FILE_alnLog") == 0 or die "aln command failed: $!";	# only for old illumina
}

#---------------------------------------
# samse
#---------------------------------------
system("bwa samse -n 100 $FILE_reference $FILE_sai $FILE_fastq > $FILE_sam 2> $FILE_samLog") == 0 or die "samse command failed: $!";


&status("convert to bam");
system("samtools view -bS $FILE_sam > $FILE_bam") == 0 or die "convert bam failed: $!";

&status("sort bam file");
my $tmp = $FILE_fastq . '.tmp';
my $tmpBam = $tmp . '.bam';
system("samtools sort $FILE_bam $tmp") == 0 or die "convert bam failed: $!";
rename $tmpBam, $FILE_bam or die "cannot find file $tmp: $!";

&status("make index");
system("samtools index $FILE_bam") == 0 or die "indexing failed : $!";






sub status{
	my ($message) = @_;
	my $time = &time();
	print "[$time] $message ... \n";
}

sub time{
	my ($sec, $min, $hour, $mday, $mon, $year) = localtime(time);
	my $time = sprintf ("%d/%0d %02d:%02d",$mon+1, $mday, $hour, $min);
	return $time;
}

